#ifndef FF_NN2_SSE_COR_TH_N_H
#define FF_NN2_SSE_COR_TH_N_H

#include <math.h>
#include <xmmintrin.h>
#include <iostream>


// my includes
#include "utils.02.h"
#include "fjf_nn2_sse_cor_n.01.h"

namespace dj01 {


class FJF_NN2_SSE_COR_TH_n {
	
	public:
	
	uint n, N;
	double *m;
	double *x0, *x1, *x2;
	double *v0, *v1, *v2;
	double *f0, *f1, *f2;
	double *j0, *j1, *j2;
	double *pot;
	double eps2;
	uint n_core;
	
	lluint nij, nij_max, nij_min; // work length for [1, nth)-threads and nth-thread
	
	dj01::FJF_NN2_SSE_COR_n forcefinder;
	
	pthread_barrier_t barrier;	
	
	struct TH {
		FJF_NN2_SSE_COR_TH_n *ff;
		double *f0,*f1,*f2, *j0,*j1,*j2, *pot;
		char *f0p,*f1p,*f2p, *j0p,*j1p,*j2p, *potp;
		pthread_t pth;
		uint id;
		uint i1_, j1_, i2_, j2_;	
		uint n1, n2;
	}*th;
	
	FJF_NN2_SSE_COR_TH_n(){
		n_core=0;
		N=n=0;
	}
		
	uint get_i02(lluint nij, uint n){
		//if (nij<1) printf("get_i()::error nij<1\n");
		llreal D = (llreal)((lluint)n*(n-1) - 2*nij) + 0.25;
		llreal k = (0.5*(2*n-1)-sqrt(D));
		return fabs(ceil(k-1));
	} 
			
	uint get_j02(uint nij, uint n, uint i){
		uint k = i+1;
		uint nijn = ((lluint)k*(k-1))/2 + (lluint)k*(n-k);
		return n-1-(nijn-nij);
	}

	uint get_i1(uint tnij_max, uint thid, uint n){
		return get_i02(tnij_max*thid+1, n);
	}

	uint get_j1(uint tnij_max, uint thid, uint n, uint i){
		return get_j02(tnij_max*thid+1, n, i);
	}
			
	static void set_ncoreopt_nijmaxmin(luint nij, uint &n_core, lluint &nij_max, lluint &nij_min){
			
		//lluint nij = ((lluint)n*(n-1))/2;
			
		lluint n_cc_min = nij/n_core;
		if (nij%n_core!=0) n_cc_min++;

		uint n_core_opt = nij/n_cc_min;		
		if (nij%n_cc_min!=0) n_core_opt++;
			
		n_core = n_core_opt;
		nij_max = n_cc_min;
		n_core>1 ? nij_min = nij - n_cc_min*(n_core_opt-1) : nij_min = 0;
	
	}
	
	// set work coordinates
	void set_work_cor(uint thid, uint n_core, uint n, uint N, lluint tnij_max, uint &i1, uint &j1, uint &i2, uint &j2){
			
		if (thid==0){
			i1 = 0;
			j1 = 1;
		} else {
			i1=get_i02(tnij_max*thid+1, N);
			j1=get_j02(tnij_max*thid+1, N, i1);
		}
		
		if (thid<n_core-1){
			i2=get_i02(tnij_max*(thid+1), N);
			j2=get_j02(tnij_max*(thid+1), N, i2);		
		} else {
			i2 = n-1; // or n-2 ?
			j2 = N-1;
		}

	}
		
	// corrects work coordinates: we need j1 to be even
	void correct_work_cor(uint N, uint &i1, uint &j1, uint &i2, uint &j2) {	// makes j1 even 
			
		if ((j2!=N-1)&& (j2%2==0)) j2++;
				
		if (j1%2!=0) 
			if (j1!=N-1) {
				if (j1!=i1+1) j1++;
			} else {
				j1=i1+2;
				i1++;
			}
	}
	
	// set sum coordinates
	void set_sum_cor(uint id, uint n_core, uint n, uint &i1, uint &i2) {
		
		if (n<n_core) n_core = n; //  reduce paralellization if n<n_core 
		
		if (id>=n_core) { i1=1; i2=0; return; } // do nothing if n>n_core and id > n
	
		uint di = ceil((double)n/n_core); // max work length
		
		i1 = id*di;
		
		if (i1>=n) { i1=1; i2=0; return; } 
			
		n-i1>di ? i2=i1+di-1 : i2 = n-1; // please think by yourself :)
		
		
	}
	
			
	static void* thread_luncher(void* any){
		((TH*)any)->ff->thread(((TH*)any)->id);
	}
	
	virtual void thread(uint id){
		
		TH thi = th[id];
				
		set_work_cor(id, n_core, n, N, nij_max, thi.i1_,thi.j1_,thi.i2_,thi.j2_);
		
		correct_work_cor(N, thi.i1_,thi.j1_,thi.i2_,thi.j2_);
		
		forcefinder.ff(n, N, 
							thi.i1_,thi.j1_,thi.i2_,thi.j2_, 
							m, 
							x0,x1,x2,
							v0,v1,v2,
							thi.f0,thi.f1,thi.f2,
							thi.j0,thi.j1,thi.j2,
							thi.pot, eps2);
		
		pthread_barrier_wait(&barrier);
		
		// combine forces, 
		set_sum_cor(id, n_core, n, thi.n1, thi.n2);
		
		double f0i,f1i,f2i, j0i,j1i,j2i, poti;
		for (uint i=thi.n1; i<=thi.n2; i++)	{
			
			f0i=f1i=f2i=j0i=j1i=j2i=poti=0;	
			for (uint j=0; j<n_core; j++) {	
				
				f0i+=th[j].f0[i]; 
				f1i+=th[j].f1[i];	
				f2i+=th[j].f2[i];
				
				j0i+=th[j].j0[i]; 
				j1i+=th[j].j1[i];	
				j2i+=th[j].j2[i];
				
				poti+=th[j].pot[i]; 
				
			}
			
			f0[i]=f0i;
			f1[i]=f1i;
			f2[i]=f2i;
			
			j0[i]=j0i;
			j1[i]=j1i;
			j2[i]=j2i;
			
			pot[i]=poti; 
		}		
		
	}
	
	virtual void ff(uint n, uint N, double *m, 
						double*x0,double*x1,double*x2,
						double*v0,double*v1,double*v2,  
						double*f0,double*f1,double*f2, 
						double*j0,double*j1,double*j2, 
						double*pot, double eps2, uint n_core) {
		
		
		if (this->N!=N || this->n_core!=n_core)   // recreate x0s.., th[i].f0.. if N or n_core changed
			if (this->n_core)
				for (uint i=0; i<this->n_core; i++) {
					
					delete [] th[i].f0p;
					delete [] th[i].f1p;
					delete [] th[i].f2p;
					 
					delete [] th[i].j0p;
					delete [] th[i].j1p;
					delete [] th[i].j2p; 
					
					delete [] th[i].potp; // delete old array
				}

		if (this->n_core!=n_core)	{
			if (this->n_core) delete [] th;
			th = new TH[n_core];
			
			for (uint i=0; i<n_core; i++) {
				th[i].id = i;
				th[i].ff = this;
			}
			pthread_barrier_init(&barrier, NULL, n_core);
		}

		if (this->N!=N || this->n_core!=n_core) 
			for (uint i=0; i<n_core; i++) {
			
				for (uint i=0; i<=4; i++) {
					char *testp; 
					double *test  = (double*)new_align(N*14, 8, 16, testp);
					delete [] testp;
				}
			
			
				th[i].f0 = (double*)new_align(N+2, 8, 16, th[i].f0p); //*10 -performance reason
				th[i].f1 = (double*)new_align(N+2, 8, 16, th[i].f1p);
				th[i].f2 = (double*)new_align(N+2, 8, 16, th[i].f2p);
				
				th[i].j0 = (double*)new_align(N+2, 8, 16, th[i].j0p); //*10 -performance reason
				th[i].j1 = (double*)new_align(N+2, 8, 16, th[i].j1p);
				th[i].j2 = (double*)new_align(N+2, 8, 16, th[i].j2p);
				
				th[i].pot = (double*)new_align(N+2, 8, 16, th[i].potp);
			}
		
		
		this->n_core = n_core;
 		this->N = N;
		this->n = n;
		this->m = m;
		
		this->x0 = x0;
		this->x1 = x1;
		this->x2 = x2;
		
		this->v0 = v0;
		this->v1 = v1;
		this->v2 = v2;
		
		this->f0 = f0;
		this->f1 = f1;
		this->f2 = f2;
		
		this->j0 = j0;
		this->j1 = j1;
		this->j2 = j2;
		
		this->pot = pot;
		
		this->eps2 = eps2;
		
		nij = (luint)N*n-((lluint)n*(n+1))/2;
		set_ncoreopt_nijmaxmin(nij, n_core, nij_max, nij_min);
		
		for (uint i=0; i<n_core; i++)
			pthread_create(&th[i].pth, NULL, &thread_luncher, (void*)&th[i]);
	
		for (uint i=0; i<n_core; i++)
			pthread_join(th[i].pth, NULL);

	}
};

}
#endif
